We  describe  a  procedure  for  the  determination  of  the  roots  of  functions  satisfying  second- 
order  ordinary  differential  equations,  including  the  classical  special  functions.  The  scheme 
is  based  on  a  combination  of  the  Priifer  transform  with  the  classical  Taylor  series  method  for 
the  solution  of  ordinary  differential  equations,  and  requires  0(1)  operations  for  the  deter¬ 
mination  of  each  root.  When  the  functions  in  question  are  classical  orthogonal  polynomials 
(Legendre,  Hermite,  etc.),  the  techniques  presented  here  also  provide  tools  for  the  evalua¬ 
tion  of  the  weights  for  the  associated  Gaussian  quadratures.  The  performance  of  the  scheme 
for  several  classical  special  functions  (prolate  spheroidal  wave  functions,  Bessel  functions, 
and  Legendre,  Hermite,  and  Laguerre  polynomials)  is  illustrated  with  numerical  examples. 


A  Fast  Algorithm  for  the  Calculation  of  the  Roots  of 
Special  Functions 


Andreas  Glaser,  Xiangtao  Liu,  and  Vladimir  Rokhlin 
Technical  Report  YALEU/DCS/TR-1367 
August  23,  2006 


The  authors  were  supported  in  part  by  the  U.S.  Department  of  Defense  under  DARPA 
Grant  #FA9550-06-l-0239  and  AFOSR  Grant  #FA9550-06-l-0197. 

Approved  for  public  release:  Distribution  is  unlimited. 

Keywords:  Roots,  Special  Functions,  Gaussian  Quadratures,  Sturm- Liouville 


Report  Documentation  Page 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 


1.  REPORT  DATE 

23  AUG  2006 


2.  REPORT  TYPE 


3.  DATES  COVERED 

00-00-2006  to  00-00-2006 


4.  TITLE  AND  SUBTITLE 

A  Fast  Algorithm  for  the  Calculation  of  the  Roots  of  Special  Functions 


6.  AUTHOR(S) 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Yale  University, Department  of  Computer  Science, New  Ha ven,CT, 06520 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


5a.  CONTRACT  NUMBER 


5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

5d.  PROJECT  NUMBER 


5e.  TASK  NUMBER 


5f.  WORK  UNIT  NUMBER 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 


13.  SUPPLEMENTARY  NOTES 


14.  ABSTRACT 


15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF: 


a.  REPORT 

unclassified 


b.  ABSTRACT 

unclassified 


c.  THIS  PAGE 

unclassified 


17.  LIMITATION  OF 

18.  NUMBER 

ABSTRACT 

OF  PAGES 

Same  as 

20 

Report  (SAR) 

19a.  NAME  OF 
RESPONSIBLE  PERSON 


Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


1  Introduction 


Zeroes  of  classical  special  functions  play  an  important  role  in  the  mathematical 
sciences,  being  related  to  Gaussian  quadratures  (both  classical  and  generalized), 
resonances  in  mechanical  and  electrical  systems,  quantum  mechanical  calculations, 
etc.  As  a  result,  algorithms  for  the  determination  of  roots  of  functions  of  this  type 
are  quite  well  developed;  most  of  them  are  based  on  the  existence  of  three-term  re¬ 
currence  relations  for  the  classical  special  functions,  and  the  associated  connection  to 
the  eigenvalues  of  certain  tridiagonal  matrices  (see,  for  example,  [6,8,2]).  The  latter 
are  dealt  with  via  well-understood  and  reliable  techniques,  such  as  the  QR  iteration. 

While  such  methods  are  quite  efficient  for  small-scale  problems,  they  require 
order  n2  operations,  where  n  is  the  the  number  of  roots  to  be  computed.  Relatively 
recently,  algorithms  for  the  diagonalization  of  tridiagonal  matrices  were  introduced 
with  CPU  time  requirements  proportional  to  n-log(n)  (see  [7]);  a  detailed  description 
of  the  resulting  fast  root  finding  scheme  is  given  in  [12].  In  this  paper  we  present 
a  scheme  which  is  unrelated  to  the  three  term  recurrence  relations  and  tends  to  be 
significantly  faster  for  large  n. 

The  subject  of  this  paper  is  the  fact  that  the  roots  of  classical  special  functions 
can  be  found  for  a  cost  proportional  to  n,  using  the  fact  that  such  functions  satisfy 
ordinary  differential  equations  of  the  form 

P(x)  ^0*0  +  q{x)  ^( x )  +  r(x )  u(x)  =  0,  (1) 

where  p(x),  q(x)  and  r(x)  are  polynomials  of  degree  2.  Using  classical  tools  such  as 
Taylor  series  and  the  Priifer  transformation,  our  approach  leads  to  remarkably  simple 
and  efficient  schemes  for  the  construction  of  classical  Gaussian  quadratures,  roots 
of  Bessel  and  prolate  spheroidal  wave  functions,  and  the  roots  of  other  functions 
satisfying  a  second-order  differential  equation  with  polynomial  coefficients. 


2  Mathematical  and  numerical  preliminaries 

In  this  section  we  summarize  several  well-known  facts  to  be  used  in  this  paper;  all 
of  these  can  be  found  in  [1,  3,  4,  11,  5].  Throughout  this  paper,  the  derivative  of  a 
function  /  :  R  — ►  R  will  be  denoted  by  f  and  p,  q,  r  will  denote  the  second-order 
polynomials  introduced  in  equation  (1).  All  functions  are  assumed  to  be  sufficiently 
smooth. 
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2.1  The  Priifer  transform 


Given  equation  (1)  and  a  differentiable  positive  function  7  :  E  — »  R,  we  define  the 
function  9  :  R  — >  R  as  the  solution  of  the  differential  equation 


9' 


sin2(0)  —  —  cos2(0) 
V  1 


Y  q  —  p'\  sin(20) 
7  P 


(2) 


where  9,^,p,q,r  are  functions  of  integrating  both  sides  of  (2)  with  respect  to  x 
and  carrying  out  elementary  manipulations,  we  obtain 


9{x)  —  arctan 


/  1  p(x)u'(x)\ 
\7(*)  «(*)  /  ’ 


(3) 


where  u  is  the  function  defined  by  equation  (1);  we  observe  that,  for  any  x  such  that 
u(x )  =  0, 

d{x)  =  (n+  1/2)tt,  (4) 

with  n  an  integer.  Similarly,  for  any  x  such  that  u'(x)=0, 


where  again  n  is  an  integer. 


9{x)  =  me, 


(5) 


Observation  1  It  often  happens  that  the  function  9  is  better  behaved  (both  nu¬ 
merically  and  analytically)  than  the  function  u  that  gave  rise  to  it.  In  such  cases, 
it  becomes  attractive  to  solve  equations  (4)  and  (5)  instead  of  equations 


u(  x)  =  0, 

(6) 

u'(x)  =  0, 

(7) 

respectively.  This  observation  is  not  new,  and  in  fact  the  above  construction  is  a 
variant  of  the  classical  Priifer  transform,  which  can  be  found,  for  example,  in  [3]. 

The  choice  of  the  function  7  in  (3)  influences  the  behavior  of  the  function  9, 
without  changing  the  solutions  of  equations  (4),  (5)  (the  latter  are  also  the  solutions 
of  equations  (6),  (7),  respectively,  with  u  being  independent  of  7).  For  example,  if 
7=1  and  equation  (1)  is  self-adjoint  (i.e.,  q  =  p') ,  then 

9'  =  —  -  sin2($)  —  r  cos2  (9)  (8) 

p 

and  the  function  9  is  strictly  decreasing  if  p  and  r  are  everywhere  positive. 

However,  it  often  happens  that  1/p  and  r  are  of  different  orders  of  magnitude; 
this  causes  the  value  of  9'  to  vary  dramatically  on  the  interval  9  =  [ — 7t/2,  7t/2]  , 
making  it  inconvenient  when  used  as  a  numerical  tool.  Figure  1(a)  illustrates  the 
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last  point.  It  shows  0  for  a  Legendre  polynomial  of  order  Z  =  100,  with  p  =  1  —  x2, 
q  =  p'  and  r  =  l (l  +  1),  in  the  interval  between  the  polynomial’s  greatest  negative 
root  and  its  least  positive  root. 

Such  step-like  behavior  complicates  many  numerical  operations  involving  (2), 
and  can  be  controlled  via  an  appropriate  choice  of  7.  Indeed,  choosing 


7  =  y/rp 

converts  (3)  into 

1  p{x)u'(x)  \ 

y/r(x)p(x)  u(x)  J 

and  substituting  (9)  into  (2)  yields 


0(x)  =  arctan 


[7  r'p  —  p'r  +  2  rq 
y  p  2  rp 


sin(2$) 

2 


(9) 

(10) 


(11) 


In  this  case,  the  large  values  of  the  product  rp  do  not  cause  excessive  variations  in 
the  magnitude  of  O',  which  is  illustrated  by  Figure  1(b),  which  shows  0  for  the  same 
Legendre  polynomial  as  in  Figure  1(a),  but  with  7  defined  in  (9). 

Whenever  4^  <  0,  the  ODE  (11)  can  be  rewritten  in  the  form 


dx  _  /  fr  r'p  —  p'r  +  2rq  sin  (20)  \  1 

dO  V  V  p  2  rp  2  / 


where  0, 7,  r,  q,p  are  functions  of  x,  and  equation  (12)  can  be  viewed  as  a  differential 
equation  defining  x  as  a  function  of  6.  Clearly,  O'  is  guaranteed  to  be  negative 
whenever 


r'p  —  p'r  +  2  rq 
Arp 


(13) 


and  the  latter  condition  is  easy  to  test  for,  and  is  satisfied  by  all  standard  special 
functions. 


2.2  Taylor  series  for  special  functions 


As  is  well  known,  any  function  u  :  R  — »  R,  that  is  sufficiently  smooth  in  the 
neighborhood  of  the  point  x0  £  R  can  be  approximated  by  its  Taylor  series 

u(x0  +  h)  =  ^  U  hk  +  e’  (14) 

k=  0  K" 

with 


|e|<  sup  (ym+1^.hw+1| , 

\x-xo\<h  l  +  !)•  / 


(15) 
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where  denotes  the  kth  derivative  of  u\  the  Taylor  series  for  v!  is 


u 


k=  1 


with 


|e|  <  sup 

\x—xo\<h  V. 


um+1(x ) 
ml 


hn 


(16) 


(17) 


If  u  satisfies  equation  (1),  and  p,  q,  r  are  polynomials  of  order  2,  then  the  con¬ 
secutive  derivatives  of  u  can  be  obtained  via  a  simple  recursion.  Indeed,  differenti¬ 
ating  (1)  A;  times  and  carrying  out  elementary  manipulations,  we  have 


pu 


(k+ 2) 


(kp1  +  q)  u(fc+1)  -  [~—^p”  +  kq'  +  r'j  u 


k(k  ^)„//„,(fc- 


(18) 


for  any  k  >  0. 


2.3  Numerical  tools 

In  this  section  we  summarize  several  numerical  techniques  to  be  used  in  the  paper. 
All  of  them  can  be  found,  for  example,  in  [4]. 

2.3.1  A  second-order  Runge-Kutta  method 

The  second-order  Runge-Kutta  method  (also  known  as  the  Heun  or  midpoint  method) 
solves  the  initial-value  problem 


y(x  0)  =  y0, 

y\x)  =  f(x,y),  (19) 

on  the  interval  [x0,x0  +  L]  by  taking  a  sequence  of  n  steps: 

Xj -f~ i  X-i  4-  h, 

kx  =  hf(xi,yi), 

k2  =  hf(xi  +  h,y  +  ki), 

y%+i  —  y%  + -(k\  +  k2) .  (20) 

where  h  =  L/n.  The  cost  of  this  algorithm  is  of  order  n,  and  the  global  truncation 

error  is  of  order  h2. 
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2.3.2  The  Taylor  series  method  for  the  solution  of  ODEs 

The  Taylor  series  method  is  a  classical  numerical  scheme  for  the  solution  of  ODEs. 
Given  equation  (1),  the  values  of  u(xk),u'(xk )  at  the  point  xk,  and  an  appropri¬ 
ately  chosen  integer  m  >  2,  the  Taylor  series  method  calls  for  repeated  differentia¬ 
tion  of  (1),  resulting  in  the  values  of  u'"{xk),u{/i) (xk), . . .  The 

latter  are  used  to  evaluate  an  approximations  to  u(xk+i),u'(xk+ 1)  via  the  expan¬ 
sions  (14),  (16),  with  h  —  Xk+i  —  Xk-  Obviously,  the  truncation  error  of  one  step 
of  the  Taylor  series  method  is  hm,  and  its  cost  is  determined  by  the  cost  of  finding 
the  derivatives  u'"(xk), . . . , u{m~^{xk),u^m){xk)]  in  the  case  of  (1),  the  recursion  (18) 
permits  u"'(xk), . . .  ,u^m~l\xk),u^m\xk)  to  be  evaluated  in  0{m)  operations. 

By  choosing  the  order  m  of  the  expansions  (14),  (16),  one  can  control  the  conver¬ 
gence  rate  of  the  Taylor  series  method,  and  often  it  makes  sense  to  choose  very  high 
orders,  obtaining  essentially  machine  precision.  In  general,  Taylor  series  algorithms 
are  not  immune  to  stability  problems;  however,  due  to  the  following  observation  this 
issue  does  not  arise  in  the  algorithm  of  this  paper. 

Observation  2  Somewhat  surprisingly,  stability  issues  are  obviated  when  all  of  the 
points  X\,X2,...  are  roots  of  u  or  u' .  Indeed,  given  the  recursion  formula  (18),  the 
expressions  (14),  (16)  can  be  viewed  as  a  linear  mapping  M  :  M2  — >  JR2  converting 
the  pair  u(xk),u'(xk)  into  the  pair  u(xk+i),u'(xk+i );  the  stability  of  the  resulting 
ODE  solver  is  determined  by  the  eigenvalues  of  M.  However,  if 

u(xk)  =  0,  (21) 

u(xk+i)  =  0,  (22) 

then  the  mapping  M  :  R2  ->  R2  is  replaced  with  a  mapping  M  :  M1  — >  R1,  con¬ 
verting  u'(xk )  into  u'(xk+ 1).  In  other  words,  the  2-dimensional  mapping  M  has 
been  replaced  with  a  simple  scaling,  not  affecting  the  roots  of  either  u  or  of  u'.  An 
analogous  argument  applies  when  equation  (21)  is  replaced  by 

v!{xk)  =  0,  (23) 

or  equation  (22)  is  replaced  by 

u'(xk+i)  =  0,  (24) 

or  when  both  equations  are  replaced. 


2.3.3  Newton’s  method 


Given  an  initial  approximation  x0,  Newton’s  method  solves  the  equation  f(x)  =  0 
iteratively,  with  the  nth  iteration  defined  by  the  formula 


%n+l  — 


fM 

I'M' 


(25) 
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If  x  is  a  simple  root  of  /  and  x0  has  been  chosen  sufficiently  close  to  x,  then  the 
convergence  is  quadratic,  i.e.,  the  number  of  correct  digits  of  the  approximation  xn 
will  double  after  each  iteration. 


2.4  Orthogonal  polynomials,  Bessel  functions,  and  prolate 
spheroidal  wave  functions 

This  section  summarizes  certain  well-known  properties  of  several  classes  of  special 
functions,  used  in  Section  4  to  test  the  performance  of  the  algorithm  of  this  paper; 
all  of  the  facts  in  this  section  can  be  found  in  [1,  11,  5]. 

2.4.1  Orthogonal  polynomials 

Given  a  (possibly  infinite)  interval  (a,  b)  and  a  positive  function  w,  such  that  the 
product  of  w  with  any  polynomial  is  integrable  on  (a,  i>),  we  say  that  po,Pi, . . .  ,pn 
is  an  orthogonal  family  of  polynomials  on  (a,  b)  with  weight  w,  if  for  each  k  = 
1, 2, 3, ... ,  the  polynomial  pk  is  of  degree  k  and 

Pi  (x)  Pj  (x)  w(x)  dx  =  {  ®  ,  (26) 

where  a,  >  0. 

In  this  subsection  we  briefly  discuss  properties  of  Legendre,  Hermite,  and  La- 
guerre  polynomials,  which  are  classical  examples  of  orthogonal  families  of  polyno¬ 
mials.  Each  of  these  classical  families  satisfies  a  homogenous  second-order  ordinary 
differential  equation. 

The  Legendre  polynomials  PQ,  Pi, . . . ,  Pn  are  the  orthogonal  family  with  respect 
to  the  weight  function  w(x)  =  1  on  the  interval  (—1,1)  for  k  =  1, . . .  ,n.  The  nth 
degree  polynomial  Pn  of  this  family  satisfies  the  differential  equation 

(1  -  x2)P”(x)  -  2x  P'n{x)  +  n(n  +  1  )Pn(x)  =  0.  (27) 

The  Hermite  polynomials  H\,  H%, . . . ,  Hn  are  the  orthogonal  family  with  respect  to 
the  weight  function  w(x)  =  exp(—x2)  on  the  interval  (—00,00),  and  the  nth  degree 
Hermite  polynomial  Hn  satisfies  the  equation 

H"(x)  -  2x  H'n{x)  -I-  2 n  Hn(x)  =  0.  (28) 

Finally,  the  Laguerre  polynomials  Li,  L2, . . .  ,Ln  are  orthogonal  with  respect  to  the 
weight  function  w(x)  =  exp(— x)  on  the  interval  (0, 00),  and  the  nth  degree  Laguerre 
polynomial  Ln  satisfies  the  equation 

xL"(x)  +  (1  -  x)L'n{x)  +  nLn(x)  =  0.  (29) 


6 


Every  class  of  orthogonal  polynomials  satisfies  a  three-term  recurrence  relation, 
i.e.,  there  exist  real  sequences  Co,  Ci,  c2,  C3, . . . ,  do,  d\,d2,  d3, . . . ,  such  that,  for  any 
to  >  1, 

CmPm+lix)  =  (x-  dm)pm(x )  -  Cm-i  pm-i  (rr).  (30) 

Thus,  given  x  G  K,  an  integer  n  >  2,  and  the  values  po(x),pi(x),  it  follows  from  (30) 
that  the  values  p2(x),p3(x),  Pa{x),  ■  ■  . , pn (x)  can  be  evaluated  in  0(n)  operations. 
Differentiating  (30),  we  obtain  the  recurrence 


CmPm+ 1(«)  =  (*  -  dm)p'm{x)  -  Cm-\lfm-X{x)  +  Pm(x),  (31) 

which  can  be  used  for  the  efficient  evaluation  of  the  derivatives  of  p2,p3, . . .. 

Specifically,  for  the  Legendre  polynomials  the  recurrence  relations  for  Pn  and  P'n 
are  (see  22.7.10,  [1]) 

O77  _L  1  rr) 

Pn+ 1  (x)  =  r~YxPn (x)  -  (32) 

O77  -i-1  T) 

p»+.w  =  (*KV)  +  P.W)  -  (33) 

with  the  initial  values  P~\{x)  =  0,  Pq{x)  =  1,  PL\{x)  =  0,  Pq(x)  =  0;  for  the 
Hermite  polynomials  the  recurrence  relations  for  Hn  and  H'n  are  (see  22.7.13,  [1]) 


Pn+ 1  (rr)  =  2rr  Hn  (x)  -  2 n  Hn_  1  (rr) ,  (34) 

H'n+ i(rr)  =  2  (rr  H'n( rr)  +  Hn(x))  -  2 n H'n_ ,{x),  (35) 

with  the  initial  values  //_  1  (rr)  =  0,  H0(x )  =  1,  H'^x)  =  0,  //o(rr)  =  0;  and  for  the 
Laguerre  polynomials  the  recurrence  relations  for  Ln  and  L'n  are  (see  22.7.12,  [1]) 


Jn+ 1 


{x) 


2n  +  1  —  rr 


n  +  1 


il+lW  = 


2n  +  1  —  rr 
n  +  1 


L'ni  rr) 


1 


(36) 

j -Ln(x)  n+lL'n- 1(»), 

(37) 

with  the  initial  values  L.^a:)  =  0,  L0( rr)  =  1,  L'_1(x)  =  0,  L(|(rr)  =  0. 

The  Legendre  and  Hermite  polynomials  of  even  degrees  are  even  and  have  a  local 
extremum  at  rr  =  0;  for  odd  degrees  they  are  odd  and  have  a  root  at  rr  =  0.  We 
will  also  make  use  of  the  fact  that  the  first  (smallest)  root  rr"  of  the  nth  Laguerre 
polynomial  satisfies 

rr"  >  — — - — — .  (38) 

Orthogonal  polynomials  are  a  classical  tool  for  the  design  of  numerical  integration 
rules  known  as  Gaussian  quadratures.  Given  the  roots  rrj , . . . ,  rrn  of  the  nth  degree 
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polynomial  pn  of  an  orthogonal  family  on  the  interval  (a,b)  with  weight  function  w, 
there  exist  real  numbers  w™, . . .  ,w%  such  that  the  formula 


[  f{x )  w(x)  dxxi^wt  f(xk ) 
Ja  k= 1 


(39) 


is  exact  for  any  /  that  is  a  polynomial  of  degree  2n  —  1  or  less. 

The  quadrature  weights  for  the  Legendre  polynomials  are  (see  15.3.1,  [11]) 


w'i 


(i-'iWW’ 

the  Hermite  quadrature  weights  are  (see  15.3.6,  [11]) 

_  tt22  n+1nl 

oii'“  —  _____________ 

‘  ’ 

and  the  quadrature  weights  for  the  Laguerre  case  are  (see  15.3.5,  [11]) 

™k 


(40) 


1 


xk  L'£{xk) ' 


(41) 


(42) 


Remark  3  The  Legendre,  Hermite  and  Laguerre  polynomials  satisfy  the  bounds 
(see,  for  example,  22.14,  [1]) 

\Pn(x)\  5;  1)  (43) 

\Hn{x)\  <^2n'2VW\ex2'\  (44) 

\Ln(x)\  <  tx'\  (45) 

The  bounds  (44)  and  (45)  indicate  that  the  Hermite  and  the  Laguerre  polynomials 
are  likely  to  cause  problems  when  used  as  a  numerical  tool.  Even  for  moderate 
values  of  x  their  values  can  become  so  large  that  it  exceeds  the  maximum  value 
representable  in  standard  floating-point  arithmetic. 

This  overflow  problem  can  easily  be  avoided  by  working  with  the  functions 

Hn{x)  =  i  Hn(x),  (46) 

7r4  2n/2vm 

Ln{x)  =  e~x/2Ln(x),  (47) 

called  normalized  Hermite  and  Laguerre  functions,  instead  of  working  with  the  poly¬ 
nomials  directly.  The  functions  Hn(x)  and  Ln( x)  have  the  same  roots  as  Hn(x)  and 
Ln(x),  and  (44),  (45)  guarantee  that  they  are  bounded  by  one.  Equations  (28),  (29) 
imply  the  following  differential  equations  for  Hn(x)  and  Ln (x) 

H"Xx)  +  (2 n  +  1  -  x2)Hn(x )  =  0,  (48) 
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x2K{x)  +  xL'n(x)  -  (lx2-(n  +  0  x^j  Ln(x)  =  0.  (49) 

The  functions  Hn  can  be  evaluated  via  the  recurrence 

Hn+ i(x)  —  Hn(x)  —  \j n  j  Hn-i{x),  (50) 

with  the  initial  values  H-\  =  0,  Hq  =  ir~^e~x2^2.  Furthermore,  if  x  is  a  root  of  Hn, 
then  H’n(x )  can  be  evaluated  via  the  recurrence 


^+l0r)  =  V  ^+T  (x"’n{x)  +  n n{x) )  "  V  ^TT^-l(x)’ 


(51) 


with  the  initial  values  H-.x  =  0,  H0  =  0.  Combining  (46)  with  (41),  we  observe  that 


Wu  = 


2e~xk 


H’2(xk) 

The  functions  Ln  can  be  evaluated  via  the  recurrence 

~  2n+  1  -x~  n  ~ 

-^n+i(#)  _  ,  t  Ln[x)  ^  (  1  Z/n_i(x), 


(52) 


n  + 1 


n+  T 


(53) 


the  same  recurrence  relation  as  for  Ln,  with  the  initial  values  L_ x(x)  =  0,  L0(x)  = 
e-x/2  If  x  is  a  root  of  Ln,  then  Ln(x)  can  be  evaluated  via  the  recurrence 


2n  +  1  -  x 
77+1 


(54) 


with  the  initial  values  L-X(x)  =  0,  L0(x )  =  0,  L'_x{x)  =  0,  L'0{x )  =  0.  Finally, 
combining  (47)  and  (42),  we  obtain 


o  %k 


(55) 


2.4.2  Bessel  functions 

The  nth  order  Bessel  function  Jn  of  the  first  kind  satisfies  the  differential  equation 

x2J'n(x)  +  x  J'n(x )  +  ( x 2  -  n2)Jn(x)  =  0.  (56) 

Figure  2  shows  a  plot  of  J50.  The  Bessel  function  of  order  n  starts  to  oscillate 
around  x  —  n,  so  that  with  xx  denoting  the  smallest  positive  root  of  Jn{x)  (see  [1], 
9.5.14) 

xx  >  n  +  rzi.  (57) 
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For  large  values  of  x,  the  asymptotic  behavior  of  Jn  is  given  by  (see  [1],  9.2.5) 

j”w“V^cK*~(n+O0+K^)'  (58) 

Remark  4  There  exist  various  classical  schemes  for  the  evaluation  of  Jn{x)  in  O(n) 
operations.  See,  for  example,  9.12,  [1]  for  Jn(x)  and  9.1.27,  [1]  for  J'n(x). 

2.4.3  Prolate  spheroidal  wave  functions 

The  prolate  spheroidal  wave  functions  (PSWFs)  corresponding  to  the  band-limit  c 
are  the  eigenvectors  of  the  integral  operator  Gc  defined  via  the  formula 

Fcmx)  =  J\icxtmdt)  (59) 

they  are  also  the  eigenvectors  of  the  differential  operator 

Gc(ip)(x)  =  (1  -  x2)V’"(x)  -  2x'tp'{x)  —  c2x2rip(x).  (60) 

The  nth  eigenvalue  of  the  operator  (60)  is  normally  denoted  by  Xn,  so  that  the  nth 
eigenvector  yjrn  of  (60)  satisfies  the  ODE 

(1  -  x2)^[x)  -  2x'ip'n(x)  +  (xcn  -  cV)i(x)  =  0,  (61) 

which  provides  the  standard  tools  for  the  evaluation  of  the  functions  ipn,  eigenval¬ 
ues  Xn>  aRd  related  quantities. 

The  function  ipn  has  n  roots;  for  even  n  the  function  ipn  is  even  and  has  a  local 
extreme  value  at  x  =  0;  for  odd  n  the  function  yjn  is  odd  and  has  a  root  at  x  =  0. 
We  refer  the  reader  to  [13,  10]  for  more  detailed  information  about  PSWFs. 

3  Algorithm 

In  this  section  we  describe  an  algorithm  for  computing  the  roots  of  a  function  u 
which  is  a  solution  of  equation  (1).  We  give  an  informal  outline  in  the  section 
below,  followed  by  a  detailed  description  in  Section  3.2. 

3.1  Informal  description 

The  algorithm  of  this  paper  finds  all  roots  of  a  function  u  satisfying  equation  (1) 
recursively  by  calculating  each  successively  bigger  root  from  the  previous  one. 
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3.1.1  First  root 


Given  a  starting  value  xs,  the  algorithm  calculates  the  smallest  root  xx  that  is  larger 
or  equal  to  xs  in  two  steps.  First,  it  finds  an  approximation  xx  to  the  initial  root  by 
calculating  9q  =  9{xs)  via  equation  (10)  and  solving  the  differential  equation  (12) 
with  the  initial  condition  x(90)  =  xs  in  the  interval  9  G  (0O,  —  7t/2)  via  the  Runge- 
Kutta  method  (20).  Due  to  (4),  the  value  of  x  at  -it/ 2  is  a  root  of  u. 

Remark  5  In  our  experience,  ten  steps  of  the  Runge-Kutta  method  are  sufficient 
to  yield  Xi  with  two  to  three  digits  of  accuracy.  Using  a  higher  order  scheme  for 
the  solution  of  the  differential  equation  (12)  and  increasing  the  number  of  nodes 
in  the  discretization  of  the  interval  (tt/2,  -tt/2),  one  could  obtain  the  root  xx  with 
arbitrarily  high  precision.  While  this  scheme  could  be  continued  to  find  all  roots 
xx,x2,  ■  ■  ■  ,xn  for  a  cost  proportional  to  n,  its  actual  cost  would  be  much  greater 
than  the  cost  of  the  scheme  of  this  paper,  where  we  solve  the  equation  (12)  to  low 
precision,  and  refine  the  solution  via  the  Newton  procedure  described  below. 

The  second  step  uses  the  Newton  procedure  (25)  to  obtain  the  root  xx ,  starting 
with  the  approximation  xx.  Since  the  Newton  process  is  quadratically  convergent, 
it  takes  about  three  steps  per  root  to  obtain  15-digit  accuracy,  and  each  of  the  steps 
requires  the  evaluation  of  u  and  vl  at  a  point  in  the  vicinity  of  xi+x. 

The  evaluation  of  u  and  v!  are  performed  by  a  scheme  specific  to  the  function  u. 
If  the  evaluation  cost  is  k  operations  then  the  cost  for  calculating  the  first  root  is 
of  order  k.  For  orthogonal  polynomials,  for  example,  u  and  u'  can  be  evaluated  via 
the  recurrence  relations  (30),  (31)  for  a  cost  which  is  of  order  of  the  degree  of  the 
polynomial. 

Remark  6  In  cases  like  the  Legendre  or  the  Hermite  polynomials,  where  an  initial 
root  or  an  extreme  value  can  be  obtained  via  symmetry  considerations,  the  compu¬ 
tation  of  the  first  root  can  be  simplified  substantially.  These  details  are  discussed 
in  Sections  3.2,  4  below. 

3.1.2  Subsequent  roots 

Given  a  root  xx  of  u,  the  algorithm  finds  the  next  bigger  root  xl+x  of  u  in  two  steps. 
Analogous  to  the  computation  of  the  first  root,  the  first  step  finds  an  approximation 
xi+x  by  solving  the  differential  equation  (12)  via  the  Runge-Kutta  scheme  (20)  on 
the  interval  9  G  (tt/2,  —tt/2)  with  the  initial  condition  x(tt/2)  =  Xi .  Due  to  (4),  xi+x 
is  an  approximation  to  the  root  Xi+X . 

The  second  step  uses  the  Newton  procedure  (25)  to  increase  the  accuracy  of 
the  approximation  xi+1.  However,  now  the  evaluation  of  u  and  vl  is  performed  via 
expansion  (14),  where  the  required  derivatives  of  u  are  computed  via  the  recur¬ 
rence  (18). 
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Remark  7  The  required  number  of  terms  in  the  Taylor  expansion  has  been  deter¬ 
mined  numerically.  For  all  examples  shown  in  this  paper,  30  terms  are  sufficient 
for  double  precision  (15-digit),  and  60  terms  are  sufficient  for  extended  precision 
(30-digit). 

The  Taylor  expansion  enables  the  evaluation  of  u  and  u'  at  a  reasonably  small 
constant  cost.  Adding  up  the  costs  for  all  n  roots,  we  observe  that  the  total  cost  of 
the  algorithm  is  of  the  order  n  +  k  operations,  where  k  is  the  cost  for  evaluating  u 
and  u'. 

The  repeated  evaluation  of  u  via  its  Taylor  expansion  can  be  seen  as  an  instance 
of  the  Taylor  series  method.  Since  the  steps  are  taken  from  one  root  to  another, 
Observation  2  obviates  the  issue  of  numerical  stability  of  the  procedure. 

Remark  8  The  algorithm  can  easily  be  adapted  to  solve  the  equation 


u'(x)  =  0  (62) 

instead  of  u(x)  =  0,  facilitating  the  computation  of  n  local  extreme  values  of  u  in 
order  n  operations. 

Remark  9  In  the  algorithm  described  above,  Taylor  expansions  are  used  to  evalu¬ 
ate  u  close  to  its  roots.  Of  course,  Taylor  expansions  can  also  be  used  to  evaluate  u 
between  roots,  thus  providing  a  scheme  for  the  fast  evaluation  of  u.  Specifically, 
given  a  set  of  n  points  X\,X2,  ■  ■  ■  ,xn,  where  Xj_i  <  xt,  and  m  roots  x\ ,  x2, . . . ,  xn  of 
u,  where  Xi-\  <  x,L,  the  scheme  finds  the  closest  root  to  each  point  X{  and  evaluates  u 
at  Xi  via  a  Taylor  expansion  at  that  root.  Including  the  cost  for  computing  the  m 
roots  of  u,  the  total  cost  of  this  scheme  is  of  the  order  of  m  +  n  operations. 

3.2  Detailed  description 

In  this  section  we  present  the  algorithm  in  some  detail.  Algorithm  1  describes  the 
core  algorithm.  As  input,  it  requires  the  smallest  root  which  is  to  be  computed  and 
the  number  N  of  roots  to  compute.  It  returns  the  roots  x\...xn  and  the  derivatives 
u'(x i) . . .  u'(x_m)  for  a  cost  of  order  N  operations.  Algorithms  2  and  3  are  auxiliary 
schemes  for  finding  the  initial  root. 

Algorithm  1 

Comment  [The  input  parameters  are  the  polynomial  coefficients  of  p,  q,  and  r  in 
equation  (1),  an  initial  root  x\  of  u,  the  number  of  roots  N  which  are  greater  than 
or  equal  to  xly  and  the  derivative  at  the  initial  root  u'{x i).  Algorithm  1  returns  the 
N  roots  of  u  which  are  greater  than  or  equal  to  X\  in  the  vector  roots.  Furthermore, 
it  returns  the  value  of  u'  at  the  roots  in  the  vector  ders. ] 
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Set  roots(  1)  =  X\  and  ders(  1)  =  u'(x\). 
do  i  =  1 ,  N  —  1 

1.  Use  the  Runge-Kutta  method  (20)  to  solve  equation  (12)  for  0  =  7r/2  to  —  n/2 
with  initial  conditions  60  =  n/2  and  x(0o)  =  roots (i).  The  resulting  value  for 
x(—n/2)  is  the  first  approximation  to  roots(i  +  1). 

2.  Use  recurrence  relation  (18)  to  calculate  the  first  30  derivatives  of  u  at  x  = 
roots(i).  Start  the  recurrence  with  u(x)  —  0,  u'{x)  =  ders(i). 

3.  Increase  the  precision  of  roots(i  4-  1)  =  x(—n/2)  via  Newton’s  method  (25), 
where  u  and  u'  are  evaluated  via  the  first  30  terms  of  the  Taylor  expansion 
around  roots(i).  Set  ders(i  +  1)  =  u'(roots{i  +  1)). 

end  do 

If  a  local  extremum  of  u  is  given  an  analogous  procedure  to  steps  1-3  of  Algo¬ 
rithm  1  can  be  used  to  calculate  the  next  biggest  root.  This  procedure  is  described 
in  Algorithm  2. 

Algorithm  2 

Comment  [The  input  parameters  are  the  polynomial  coefficients  of  p,  q,  and  r  in 
equation  (1),  an  extreme  value  xe  of  u,  and  u(xe).  Algorithm  2  returns  X\,  the 
smallest  root  of  u  which  is  bigger  than  xe,  and  the  derivative  v!(x\)] 

1.  Use  the  Runge-Kutta  method  (20)  to  solve  equation  (12)  for  9  =  0  to  —n/2 
with  initial  conditions  90  =  0  and  x(90)  =  xe.  The  resulting  value  for  x(-n/2) 
is  the  first  approximation  to  x\. 

2.  Use  recurrence  relation  (18)  to  calculate  the  first  30  derivatives  of  u  at  xe. 
Start  the  recurrence  with  u(xe),  u'(xt)  —  0. 

3.  Increase  the  precision  of  X\  =  x(—n/2)  via  Newton’s  method  (25),  where  u 
and  v!  are  evaluated  via  the  first  30  terms  of  the  Taylor  expansion  around  xe. 

If  the  first  roots  of  u  lie  in  a  regime  where  the  Taylor  expansion  has  a  poor 
convergence  rate  (see,  for  example,  the  results  on  Laguerre  polynomials  in  Section 
4.3),  the  following  algorithm  can  be  used. 

Algorithm  3 

Comment  [Algorithm  3  requires  a  routine  for  evaluating  u  and  u' .  Its  input  pa¬ 
rameters  are  the  polynomial  coefficients  of  p ,  q,  and  r  in  equation  (1),  and  a  starting 
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value  xs  which  is  less  than  the  smallest  root.  Algorithm  3  returns  x1;  the  smallest 
root  of  u  which  is  bigger  than  xs  and  u'(x i).] 

1.  Compute  9(xs )  via  equation  (10). 

2.  Use  the  Runge-Kutta  method  (20)  to  solve  equation  (12)  for  9  =  9{xs)  to 
— 7t/2  with  initial  conditions  9q  =  9{xs)  and  x(90)  =  xs.  The  resulting  value 
for  x(—tx/2)  is  the  first  approximation  to  X\. 

3.  Increase  the  precision  of  x\  =  x(—tx/2)  via  Newton  iterations,  where  u  and  u' 
are  evaluated  by  the  given  routine,  which  is  specific  to  u.  Return  x\  and  u'(x\). 


4  Implementation  and  numerical  results 


We  have  implemented  the  algorithm  of  this  paper  in  FORTRAN  77  and  tested  it  on 
several  classical  families  of  special  functions.  Below  we  discuss  our  implementation 
for  Legendre,  Hermite,  and  Laguerre  polynomials,  Bessel  functions,  and  prolate 
functions.  For  each  example  we  give  timings  and  accuracy  for  different  problem 
sizes.  The  computation  time  has  been  measured  on  an  Intel  Pentium  4,  2.4  GHz 
with  1  GB  RAM.  The  given  accuracy  is  based  on  computations  in  double  precision. 

The  accuracy  shown  in  the  tables  has  been  measured  by  computing  the  roots 
in  double  (15-digit)  and  in  extended  (30-digit)  precision.  Given  N  roots  x[d> . .  .xff 
of  the  function  u  in  double  precision  and  the  exact  values  (calculated  in  extended 
precision)  x[9^ , . . . ,  x$ ,  the  absolute  accuracy  ra  and  the  relative  accuracy  rr  are 
computed  as 


ra 


max 

i€{  1 . N} 


(63) 


ry 


max 
<€{1,. ...JV} 


X, 


(d)  _  X(q) 


X, 


(?) 


(64) 


Remark  10  One  of  the  principal  applications  of  the  roots  of  orthogonal  polynomi¬ 
als  is  numerical  integration,  as  given  by  formula  (39).  In  the  tables  for  the  orthogonal 
polynomials  (Tables  1,  2,  and  3),  we  show  the  relative  accuracy  of  the  roots  and  the 
absolute  accuracy  of  the  weights,  as  computed  by  the  algorithm  of  this  paper.  In 
addition,  we  provide  the  accuracy  of  the  quadrature  formula  (39). 

For  an  orthogonal  polynomial  of  degree  n  we  measure  this  accuracy  by  using  (39) 
to  integrate  the  polynomial  of  degree  \  n  of  the  same  family,  e.g.,  the  accuracy  of 
1000  computed  Laguerre  nodes  is  tested  by  integrating  the  Laguerre  polynomial 
of  degree  1500.  The  analytical  value  of  this  integral  is  zero,  so  that  the  resulting 
numerical  value  is  a  direct  measure  for  the  error. 

The  orthogonal  polynomial  which  is  to  be  integrated  is  evaluated  via  the  proce¬ 
dure  described  in  Remark  9.  To  avoid  additional  numerical  errors,  these  computa¬ 
tions  are  performed  in  extended  precision  (30  digits). 
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4.1  Legendre  polynomials 

Equation  (27)  is  the  ODE  defining  the  nth  degree  Legendre  polynomial  Pn.  Due  to 
the  symmetry  of  Legendre  polynomials  it  is  sufficient  to  compute  only  the  positive 
roots.  If  n  is  even,  Pn  has  an  extreme  value  at  x  =  0,  and  Algorithm  2  is  used  to 
find  the  first  positive  root,  followed  by  Algorithm  1  to  find  the  nj 2—1  subsequent 
roots.  If  n  is  odd,  Pn  has  a  root  at  x  =  0,  and  Algorithm  1  is  used  to  find  the 
(n  —  l)/2  positive  roots.  The  values  for  Pn( 0)  or  P^( 0)  which  are  required  as  input 
for  Algorithm  2  or  1  are  calculated  via  the  recurrence  relations  (32)  and  (33).  Finally, 
the  Legendre  quadrature  weights  are  calculated  via  formula  (40). 

Table  1  shows  the  total  computation  time,  as  well  as  the  accuracy  of  the  roots, 
the  weights,  and  the  quadrature  formulae  (see  Remark  10).  The  total  computation 
time  includes  the  computation  of  the  quadrature  weights. 

4.2  Her  mite  polynomials 

To  avoid  numerical  overflow,  we  have  applied  our  algorithm  to  the  Hermite  functions 
Hn  instead  of  the  polynomials  (see  Remark  3).  The  differential  equation  for  the 
Hermite  functions  is  given  in  (48). 

Since  Hermite  polynomials  have  the  same  symmetry  as  Legendre  polynomials, 
the  computation  of  the  first  root  and  the  application  of  the  main  algorithm  is  com¬ 
pletely  analogous.  The  values  Hn( 0)  and  H'n (0)  can  be  calculated  via  the  recurrence 
relations  (50),  (51).  Formula  (52)  is  used  to  calculate  the  Hermite  weights  from  the 
returned  roots  and  derivatives  of  the  Hermite  functions.  Table  2  shows  computation 
time  and  accuracy  (see  Remark  10). 

4.3  Laguerre  polynomials 

As  in  the  Hermite  case,  we  have  used  the  differential  equation  for  the  Laguerre 
functions  given  in  equation  (49)  instead  of  the  ODE  for  the  Laguerre  polynomials 
(see  Remark  3).  Equation  (38)  provides  the  starting  value  xs  =  4^2,  and  recurrence 
relations  (53)  and  (54)  are  used  for  the  evaluation  of  the  function  and  its  derivative. 
The  quadrature  weights  are  calculated  via  formula  (55).  Table  3  shows  timings  and 
accuracy  (see  Remark  10). 

Remark  11  The  differential  equations  (29),  (49)  have  a  singularity  at  x  =  0,  lead¬ 
ing  to  bad  convergence  rates  of  the  Taylor  expansion  of  Ln  close  to  x  =  0.  In 
order  to  avoid  this  issue,  the  first  20  roots  are  calculated  via  repeated  application 
of  Algorithm  3;  the  subsequent  roots  are  then  calculated  via  Algorithm  1. 


15 


4.4  Bessel  functions 

The  differential  equation  for  the  nth  order  Bessel  function  Jn  of  the  first  kind  is 
given  in  equation  (56).  The  first  root  is  found  via  Algorithm  3  with  the  starting 
value  xs  =  n  +  ni  as  given  by  equation  (57).  The  values  for  Jn(xs)  and  J'n{xs) 
are  calculated  by  means  of  a  standard  scheme,  which  can,  for  example,  be  found 
in  9.12,  9.1.27,  [1]. 

We  compute  the  first  20  n  roots  of  Jn.  For  larger  roots,  Jn  will  start  to  behave 
like  a  cosine,  due  to  the  asymptotic  formula  (58),  and  the  distance  between  two 
consecutive  roots  approaches  n.  Computation  time  and  accuracy  for  the  Bessel 
roots  are  shown  in  Table  4. 

4.5  Prolate  spheroidal  wave  functions 

The  differential  equation  for  prolate  spheroidal  wave  functions  is  given  in  equa¬ 
tion  (61).  We  refer  the  reader  to  [13,  10]  for  a  procedure  to  calculate  the  parame¬ 
ter  Xn ■  The  computation  of  the  first  root  and  the  application  of  the  algorithm  to 
prolate  functions  is  completely  analogous  to  Legendre  polynomials,  since  they  have 
same  symmetry.  Table  5  shows  computation  time  and  accuracy. 
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Figure  1:  Figure  (a)  shows  9  as  defined  in  equations  (2)  and  (3)  for  a  Legendre 
polynomial  of  order  100  with  7=1.  Figure  (b)  shows  9  for  the  same  polynomial 
with  7  =  rp . 


Figure  2:  Plot  of  J50 


number  of 
nodes 

relative  acc. 
of  roots 

absolute  acc. 
of  weights 

accuracy 
of  quadrature 

computation 
time  in  sec. 

4  x  10~15 

2  x  10~16 

6  x  10~16 

0.02 

1'  S 

5  x  10"15 

1  x  10"16 

4  x  10~15 

0.15 

100000 

1  x  10~16 

2  x  10~15 

1.30 

3  x  10“14 

2  x  10~15 

12.31 

Table  1:  Timings  and  accuracy  for  Legendre  roots  and  weights  (see  Remark  10). 
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number  of 
nodes 

relative  acc. 
of  roots 

absolute  acc. 
of  weights 

accuracy 
of  quadrature 

computation 
time  in  sec. 

1000 

5  x  10~15 

2  x  10~16 

0.02 

10000 

1  x  10~14 

2  x  10~15 

0.16 

100000 

5  x  10-14 

1 

6  x  10~16 

1.54 

1000000 

9  x  10~13 

2  x  10"16 

15.36 

Table  2:  Timings  and  accuracy  for  Hermite  roots  and  weights  (see  Remark  10). 


number  of 
nodes 

relative  acc. 
of  roots 

absolute  acc. 
of  weights 

accuracy 
of  quadrature 

computation 
time  in  sec. 

1  000 

4  x  10~12 

2  x  10~13 

1  x  10~13 

0.04 

10000 

X 

i — * 

o 

5 

9  x  10~13 

2  x  10~13 

0.37 

100000 

1  x  10~u 

8  x  lO"12 

3.27 

1000000 

5  x  10~u 

31.18 

Table  3:  Timings  and  accuracy  for  Laguerre  roots  and  weights  (see  Remark  10). 


number  of 
computed  roots 

order  n  of  Bessel 
function 

relative  acc. 
of  roots 

computation 
time  in  sec. 

2  000 

HEOTiBlsflB 

20000 

200000 

5.83 

Table  4:  Timings  and  accuracy  for  the  roots  of  Bessel  functions. 


number  of 
computed  roots 

band- 
limit  c 

relative  acc. 
of  roots 

computation 
time  in  sec. 

1000 

4  x  10~15 

10000 

6  x  lO" 15 

100000 

2  x  10-14 

2.33 

4  x  10“14 

22.27 

Table  5:  Timings  and  accuracy  for  the  roots  of  prolate  spheroidal  wave  functions. 
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